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Abstract 

We present the results of a large scale molecular dynamics computer simula- 
tion in which we investigated the static and dynamic properties of a silica melt 
in the temperature range in which the viscosity of the system changes from 
O(10" 2 ) Poise to O(10 2 ) Poise. We show that even at temperatures as high as 
4000 K the structure of this system is very similar to the random tetrahedral 
network found in silica at lower temperatures. The temperature dependence 
of the concentration of the defects in this network shows an Arrhenius law. 
From the partial structure factors we calculate the neutron scattering func- 
tion and find that it agrees very well with experimental neutron scattering 
data. At low temperatures the temperature dependence of the diffusion con- 
stants D shows an Arrhenius law with activation energies which are in very 
good agreement with the experimental values. With increasing temperature 
we find that this dependence shows a cross-over to one which can be described 
well by a power-law, D oc (T — T c ) 7 . The critical temperature T c is 3330 K 
and the exponent 7 is close to 2.1. Since we find a similar cross-over in the 
viscosity we have evidence that the relaxation dynamics of the system changes 
from a flow-like motion of the particles, as described by the ideal version of 
mode-coupling theory, to a hopping like motion. We show that such a change 
of the transport mechanism is also observed in the product of the diffusion 
constant and the life time of a Si-0 bond, or the space and time dependence 
of the van Hove correlation functions. 

PACS numbers: 61.20.Lc, 61.20.Ja, 02.70.Ns, 64.70.Pf 
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I. INTRODUCTION 



In recent years the dynamics of supercooled liquids has been the focus of many investiga- 
tions and due to these efforts our understanding of this dynamics has increased significantly. 
Although very recently new interesting phenomena in this dynamics have been found, such 
as the dynamical heterogeneities |IJ or the aging dynamics of the system at very low tem- 
peratures [§], which so far are not understood, we have by now at least a fair understanding 
of the relaxation dynamics at temperatures a bit above the experimental glass transition 
temperature T g . In particular it has been demonstrated that the so-called mode-coupling 
theory (MCT) ||, a theory which relates the slowing down of the dynamics upon cooling to 
nonlinear feedback effects, is able to give a qualitative and quantitative correct description 
of glass-formers in which the interaction between particles is not too different from the one 
of a hard-sphere system, such as colloidal systems |4[], the molecular glass- former orthoter- 
phenyl and even systems in which some hydrogen bonding is present, such as glycerol ||. 
The main prediction of the theory is that there exists a critical temperature T c at which 
the relaxation dynamics of the system changes qualitatively in that an unexpectedly fast in- 
crease of the local activation energy of the transport coefficients, such as the viscosity or the 
inverse of the diffusion constant, is observed. Qualitatively such a temperature dependence 
can readily be seen in the so-called Angell plot |J7|, where the logarithm of the viscosity is 
plotted versus T g /T, in that the so-called fragile glass-formers show, at around 1.1-1.3T S , 
a sharp upward bend. Because of the similarity between this bend and the prediction of 
MCT, many investigations on the dynamics of glass-formers have been done for these types 
of systems, so that the validity of the theory could be checked. Apart from this, most fragile 
systems have glass transition temperatures which are readily accessible in an experiment, 
thus making it relatively simple to investigate the dynamics at and above T g . This is not the 
case for most strong glass-formers, since their T g is often above 500 K, such as SiC>2, B 2 03, 
or GeC>2, to name a few. Because of their high value of T g the microscopic dynamics of these 
systems is understood in much less detail than the one of fragile glass-formers. Sidebottom et 
al. H have investigated the dynamics of B 2 03 by means of photon correlation spectroscopy, 
and very recently Wischnewski et al. have used neutron scattering experiments to study 
the dynamics of Si0 2 also at high temperatures. Despite these investigations our knowledge 
on the dynamics of strong glass-formers is significantly less than the one of fragile ones. 

Valuable insight into the dynamics of supercooled liquids have also been obtained through 



computer simulations. Whereas many of the earlier studies |LT| have investigated the dy- 
namics of strong glass-formers, most of the more recent studies focused on simple liquid like 
systems, such as soft spheres or Lennard- Jones systems, since in the latter type of models 
the interaction between the particles is much simpler than in the one of strong glass-formers 
and hence the system can be probed on larger time windows. Very recently however, the 
investigation of strong glass-formers, mainly silica, by means of computer simulations has 
become more popular again. From such simulations it is hoped that insight is gained into, 
e.g., the nature of the so-called boson-peak []TT|-|TB]], a dynamical feature at around ITHz 
whose origin is still a matter of debate PJIl, or to investigate the structure and dynamics 



at different pressures and temperatures on a microscopic level |]15 | 



Despite the mentioned investigations we are however still lacking a comparable detailed 
characterization of the relaxation dynamics of strong glass-formers as it has been obtained 
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for simple liquids [p!6| -[T9|| and water | 20| . In the present work we therefore present the 
results of a molecular dynamics computer simulation of silica, the prototype of a strong 
glass-former, in which such an analysis has been done. These results will allow us to gain 
more insight into the dynamics of this system, to relate it to its structural properties and 
to compare these results with the ones obtained for fragile systems. The rest of the paper 
is organized as follows: In Sec. we give the details of the model and the simulation, in 
Sec. [TTT] we present the results, and in Sec. [TV] we summarize and discuss them. 



II. MODEL AND DETAILS OF THE SIMULATIONS 

As already mentioned in the Introduction, various aspects of the dynamics of amorphous 
silica have already been investigated in the past. In these simulations potentials with dif- 
ferent levels of accuracy have been used, most of which give a satisfactory description of 
the structural (i.e. static) properties of amorphous silica. Very recently Hemmati and Angell 
have pointed out, however, that the predictions of these various potentials regarding dynam- 
ical properties, such as the diffusion constant at low temperatures, can differ by orders of 
magnitudes [21]. For the investigation of the dynamical behavior it is therefore important 



to choose the potential with particular care. One of the most reliable potentials seems to 
be the one proposed by van Beest et al. (BKS) which was developed by using a mixture of 
ab initio calculations and classical lattice dynamics simulations f22fl . The BKS potential is 
given by 

0( r . j ) = M^L + Aij eX p(- J Bi i r ii ) - ^ . (1) 
r ij ' ij 

Here is the distance between ions i and j, and the values of the partial charges qi and the 
constants Aij, B^, and can be found in Refs. p2| , p3| . As we have already done in previous 



investigations of this model |13|J23 -P7| , we have truncated and shifted the non-Coulombic 
part of this potential at a cut-off radius of 5.5A. This modification of the potential has the 
positive effect that at normal pressure the density of the glass is around 2.3g/cm 3 , a value 
that is very close to the one of real silica glass, 2.2g/ cm 3 . In previous work we have found 
that the dynamics of strong glass formers show quite significant finite size effects ||25|| . In 
order to avoid them we used a system which is relatively large for simulation of supercooled 
liquids, i.e. 8016 ions. During the runs the size of the cubic simulation box was fixed to 
L = 48.37A. Thus the density of the system was kept constant at 2.37g/cm 3 , which leads to 
a pressure which is a bit above normal pressure (discussed in more detail below). 

The Coulombic part of the potential was evaluated by means of Ewald-sums with a 
constant aL=12.82 and by using for the Fourier part of the Ewald-sum all /c-vectors 



of magnitude less than k c = 27rv51/L. A short calculation shows that, although such a 
small value of k c is sufficient to calculate the interaction between two ions to a satisfactory 
accuracy, it would be insufficient to set the absolute error of the total Coulombic energy 
below an acceptable level. However, it can be shown [3~0|, that this small value of k c leads 



just to a constant shift of the energy scale and hence does not affect the forces. 

Using the potential given by Eq. (0) and the masses of 28.086u and 15.9994u for the silicon 
and oxygen atoms, respectively, we equilibrated the system by doing a NVT simulation 
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(temperature was kept constant by means of a stochastic collision algorithm). At each 
temperature investigated, we made sure that the length of these equilibration runs exceeded 
the typical relaxation times of the system. The temperatures considered were 6100 K, 
5200 K, 4700 K, 4300 K, 4000 K, 3760 K, 3580 K, 3400 K, 3250 K, 3100 K, 3000 K, 
2900 K and 2750 K. In order to improve the statistics of the results we averaged at each 
temperature over two independent runs. After the equilibration we started microcanonical 
runs by integrating the equations of motion with the velocity form of the Verlet algorithm. 
The time step used was 1.6fs, which is sufficiently small to guarantee a negligible drift of 
the total energy of the system at high and intermediate temperatures (6100K> T > 3100K) 
over the whole run. For the lowest temperatures (T < 3000K) the runs were so long that 
the drift in the energy could no longer be neglected. Therefore we rescaled every one million 
time steps the velocity of the particles to bring the system back to the energy at which it 
was started. This frequency should be small enough to avoid any effect on the relaxation 
dynamics. At the lowest temperature the length of the runs were 12 million time steps (and 
13 millions for the equilibration), giving a total simulation time of about 20 ns. The overall 
computational effort (equilibration plus production) for this temperature was about 13 years 
of single processor time on a CRAY-T3E (or 2.5 months on 64 processors). 



III. RESULTS 

In order to understand the dynamical behavior of the system it is useful to know some- 
thing about its static properties since we will show that the transport mechanism of the 
atoms is intimately related to the network structure of the system. Therefore we present in 
the following subsection some of the static properties before we discuss in the subsequent 
subsection the dynamics of the system. 



A. Static properties 

One of the most remarkable properties of amorphous silica is that it shows a density 
anomaly at around 1820 K [E8J. Although in a constant volume simulation such an anomaly 
can of course not be observed in the density, it can easily be seen in the temperature 
dependence of the pressure p, which is shown in Fig. [I]. We see that after a fast decrease 
of p between 6100 K and 5200 K, the pressure reaches a minimum at around 4900 K, thus 
showing that the density at constant pressure has a maximum, and then increases again 
with decreasing temperature. Since the experimental value at which this density anomaly 
occurs is 1820 K [^8| we confirm the result of the constant pressure simulation f23|, that the 



BKS-potential overestimates the temperature at which the density anomaly is observed. 

That this extremum in pressure is not accompanied by a significant change in the struc- 
tural quantity g a/3 , the partial radial distribution functions between particles of type a and 
f3 [jS!], is demonstrated in Fig. [| where we show the three g a p for three different tempera- 



tures. What can also be recognized from these figures is that the BKS-model predicts that 
even at a temperature of 6100 K the arrangement of the particles shows a structure which is 
much more similar to the one of a liquid at its triple point than to the one of a gas. It should 
be realized that this result is not in contradiction to the experimental result that a silica 



4 



melt at several thousand Kelvin does evaporate, if a free surface is present. We have found 
that if the geometry of the simulation does contain such a free surface the BKS model also 
predicts that the atoms will evaporate from the system |3^]. However, due to the presence 
of the periodic boundary no free surface is present in the simulation discussed in this work 
and thus no evaporation takes place. This situation is thus similar to a melt in the interior 
of the earth. 

We also remark that in the temperature range considered, the time scale of the relaxation 
dynamics of the system changes by about four decades (this will be discussed in the next 
subsection). Thus we find that a relatively small change in the structure of the melt is 
accompanied by a dramatic change in the dynamics, an observation which is qualitatively 
similar to the one found for simple liquids [|T6| — pLQ[1 . This fact thus gives support to the 
point of view put forward by the mode-coupling theory ||, which predicts that although 
the temperature dependence of all structural quantities of a supercooled melt is very weak, 
most dynamical quantities connected to the density fluctuations of the system will show a 
very strong temperature dependence. Below we will discuss these aspects in more detail. 

From the figure we also recognize that for all temperatures there is a well defined min- 
imum between the first and second peak. Therefore it is possible to identify for each ion 
its nearest neighbors by requiring that they are within the first nearest neighbor shell, de- 
fined by the location of the mentioned minimum. For the following analysis we used the 
values 3.64A, 2.35A, and 3.2lA for the location of this minimum in the Si-Si-, Si-O-, and 
O-O-correlation. In Fig. |^ we show the probability P a p that a ion of type a has exactly z 
nearest neighbors of type (3, for all relevant values of z. From Figs. §a and |3|b we recognize 
that even at high temperatures more than 85% of the silicon and oxygen atoms are four 
and two fold coordinated, respectively, i.e. that they are in the local environment expected 
from an ideal (corner-shared) tetrahedral network. At the lowest temperature investigated 
this fraction is higher than 99%. The remaining defects are silicon atoms that are three 
and five fold coordinated and oxygen atoms that are three and one fold coordinated, the 
latter thus forming dangling bonds. We have found that for temperatures below 3700 K the 
probability for the occurrence of such defects is described very well by an Arrhenius law, 
P a /3 = ^0/3 exp(— e Qj g/T) (dashed curves). As prefactors 7r a/3 and activation energies e a p we 
found: 7r Si0 =4.4, es;o=17130 K for z = 5, 7r Si0 =58.6, esio=31100 K for z = 3, 7r si=3.2, 
eosi=17730 K for z = 3 and 7Tosi=8.9, eosi=24760 K for z = 1. From these activation ener- 
gies we recognize that at low temperatures defects in which a silicon atom is coordinated by 
z = 5 oxygen atoms or one in which an oxygen atom is bound to z = 3 silicon atoms are by 
far the most frequent ones. In Refs. [23 , 33] it was shown that once the system has reached 
the glass transition temperature, the distributions Psio an d -Posi do not change anymore, 
i.e., that at low temperatures these distributions are given by the ones at the glass transi- 
tion temperature. By using the above given activation energies we can therefore extrapolate 
the curves to the experimental value of T g and thus estimate the density of the defects of 
real silica below T g . Using the value T ff =1450 K we predict that five fold coordinated 
silicon atoms occur with probability 3.2 • 10~ 5 and three fold coordinated oxygen atoms with 
probability 1.5 • 10~ 5 . 

Whereas the temperature dependence of the Si-0 and O-Si defects is quite simple, we 
see from Figs. |3|c and |3]d that this is not the case for the defects involving the structure 
on somewhat larger length scales, since the functional form of the various curves is not 
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easy to identify. Thus the main information that can be obtained from these figures is that 
locally the system seems to approach the geometry of a random tetrahedral network. In 
this structure the most probable configuration is that an oxygen atom has six other oxygen 
atoms as its second nearest neighbors and a silicon atom has four other silicon atoms as its 
second nearest neighbors. 

The structure at larger length scales can best be studied by means of the partial structure 
factors S a p(q) [0, since any large scale feature will show up at wave- vectors q which are 
smaller than the location of the peak between neighboring particles, i.e. in our case the Si-0 
correlation. In Fig. |]we show the three S a p(q) for the same temperatures already discussed 
in Fig. 0. The peak corresponding to the nearest neighbor Si-0 distance is around 2.8A -1 . 
At smaller wave-vectors an additional peak can be seen, the so-called first sharp diffraction 
peak (FSDP). The microscopic reason for this peak is the local chemical ordering of the ions 
in tetrahedra-like structures and the location of the peak is related to the distance between 
neighboring tetrahedra. The interesting information that can be obtained from Fig. [| is that 
the FSDP is observed already at temperatures as high as 4000 K. At first glance it might 
seem a somewhat unrealistic prediction of the BKS model that the open structure of the 
network should be present even at such high temperatures. However, if one recalls that even 
at temperatures as high as 2750 K silica is a quite viscous liquid (r/ « 2000 Poise) |34|], this 
result is no longer that surprising, since such a high value of the viscosity shows that even at 
these high temperatures the particles move quite slowly, and that therefore their motion is 
strongly hindered by their neighbors. Therefore it is not surprising that the different S a p(q) 
show that there is indeed a lot of structure in the system even at high temperatures. 

Also included in Fig. |]are curves for T = 300 K. These were obtained by using the equili- 
brated configurations at T = 2900 K as starting configurations of a cooling run in which the 
temperature was decreased linearly within one million time steps to zero Kelvin. This corre- 
sponds to a cooling rate of about 1.8-10 12 K/s. With this cooling rate the system falls out of 
equilibrium at around 2850 K and thus this temperature corresponds to about the value of 
the Active temperature of the glass. The configurations that we obtained at T = 300 K were 
annealed for another 500,000 time steps before we started to calculate the structure factors. 
From the figure we recognize that the partial structure factors for this temperature seem to 
evolve very smoothly from the (equilibrium) ones at higher temperature, which however is 
not quite correct, as we will demonstrate below. 

In order to compare our results for the structure factors with the one of real silica we 
have also calculated the neutron scattering function S n (q) from 



N 



Sn ^ = at h2 T at h2 Wexp(zq • (r fc - r,))), (2) 
^VsiOsi + 7V o»o k ,i 

where k 6 {Si, 0} are the neutron scattering cross sections, and (.) is the thermal average. 
Susman et al. [[35] report for bsi and bo the values 0.4149 • 10~ 12 cm and 0.5803 ■ 10~ 12 cm, 
respectively. Using these values and Eq. (||) we obtain for T = 300 K the S n (q) shown in 
Fig. |5|. Also included in the figure is the neutron scattering result of Price and Carpenter 



at the same temperature [36]. We recognize that qualitatively, as well as quantitatively, the 
two curves agree very well, in that not only the location, but also the height of the various 
peaks is reproduced well. This shows that the BKS potential is indeed able to reproduce this 
structural quantity. A similar good agreement was obtained by Taraskin and Elliott who 
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calculated S n (q) for Si02 within the harmonic approximation One significant difference 
between the theoretical and the experimental curve can, however, be noticed, and this is 
their behavior at small wave-vectors. Whereas an extrapolation of the theoretical curve to 
q = seems to give a value around 0.16, the experimental curve seems to extrapolate to a 
significantly smaller value. It is well known that this value of the total structure factor, 

%)=iV- 1 E(exp((zq-(r fc -r,))), (3) 

k,i 

is related to the isothermal compressibility Kt of the system via 

lim S (q) = pk B Tn T , (4) 

where p = N/V is the particle density. We therefore extrapolated our partial structure 
factors to q = and thus obtained S(q) for q — > 0. This quantity is shown in Fig. ||] as 
a function of temperature. We see that as long as we are in equilibrium, T > 2750 K, 
S(q = 0) decreases continuously with decreasing temperature. At our (computer) glass 
transition T m 2850 K, however, the structure freezes in and thus at T = 300 K S(q = 0) 
is essentially given by its value at this glass transition temperature^. Thus we see that 
it is not surprising that in Fig. |5| the theoretical curve at small wave-vectors is above the 
experimental one. 

In the inset of Fig. |S| we show the compressibility of our system as calculated from Eq. 
We see that also this quantity has a maximum at a temperature which corresponds to the 
location of the density anomaly of the system, i.e. for this model at around 4700 K. No 
experimental data for in this temperature range are known to us. However, Fraser f3~? 



reports for T = 1673 K a compressibility around 2.7 ■ 10 -5 GPa -1 . If we extrapolate the 
(equilibrium) compressibility data in Fig. |5] to T = 1673 K, we obtain a kt at around 
1.35 • 10~ 5 GPa -1 which is a factor of two smaller than the one given by Fraser. Thus we 
have evidence that the BKS model is a bit too stiff. The possibility that this discrepancy is 
due to the slightly enhanced pressure of our system should, however, not be disregarded. 



B. Dynamic properties 

Having studied the static properties of the silica melt we will discuss in this subsection 
the dynamic ones. 

One of the simplest quantities to study the dynamics of a fluid system on a microscopic 
level is the mean squared displacement (MSD) (r 2 (t)) of a tagged particle (of type a), which 
is given by 

1 Nat 

(r 2 (t)) = —Y,(\n(t)-r l (0)\ 2 ) • (5) 

iV « i=i 



We remind the reader that the glass at 300 K was produced by cooling the system from 2900 K. 
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The time dependence for this quantity is shown in Fig. |7| for all temperatures investi- 
gated. For high temperatures, curves to the left, two time regimes can be distinguished. At 
short times the motion of the particles is ballistic, i.e. rj(£) ~ ^(0) + r/t, and the MSD is 
proportional to t 2 . For long times the motion is diffusive and hence we have (r 2 (t)) oc t. 
These two regimes are also seen at low temperatures, curves to the right. However, at these 
temperatures these two regimes are separated by a third one, in which the MSD does not 
change significantly over a time span which extends over several decades in time. The mi- 
croscopic reason for this behavior is the so called cage-effect, i.e. the fact that in this time 
regime the tagged particle is temporarily trapped by the particles surrounding it and hence 
make an escape from its cage very unlikely. Only for sufficiently large time the particle 
succeeds to escape this (dynamic) cage and the MSD starts again to increase significantly. 
The dynamics of the particles related to this cage is usually called the /3-relaxation and the 
mentioned mode-coupling theory makes detailed predictions about this dynamics, as it will 
be discussed in more details below. 

What is also noticeable in the curves at low temperatures is that the transition from the 
ballistic to the cage regime is quite different from the one found in a simple liquid, see e.g. 



Ref. ||19|| , in that i) immediately after the ballistic regime the curves show a small shoulder 
at around 0.03 ps and ii) a peak is observed at around 0.2 ps. The first of these features, 
which is more pronounced in the curves for the silicon atoms, is most likely the result of the 
complex local motion of the atoms in the open tetrahedral network, such as the bending and 



stretching of the tetrahedra. The second one is, as already pointed out by Angell et al. [p8 



related to the so-called boson peak, a vibrational feature at low frequencies whose precise 



origin is currently still a matter of debate |12|-|I4" 



We also point out that the above mentioned cage-effect becomes observable already 
at temperatures as high as 3580 K. Usually this effect is associated with the glass-former 
being in a supercooled state. However, since the melting temperature T m of silica is around 
2000 K, all our simulations are above T m and hence we do not investigate the system in its 
supercooled regime at all. This shows that the melting temperature is, from the point of view 
of the cage effect, an irrelevant temperature. That this somewhat surprising result is not a 
particularity of the silica model studied here can be inferred from the fact that experimental 
studies of glycerol and B 2 03 also show a strong non-Debye relaxation at temperatures well 
above their melting temperature P,|39"f . All these results give support to the point of view of 
the earlier mentioned MCT, that the whole slowing down of the system is a purely kinetic 
phenomenon and thus not related to any thermodynamic singularity of any kind. 

Using the Einstein relation limt^ 00 (r 2 (t)/6t) = D, it is easy to calculate the diffusion 
constants D from the MSD. These are plotted in Fig. [5] as a function of the inverse temper- 
ature. As expected for a strong glass-former, we find that at low temperatures the diffusion 
constants show an Arrhenius dependence. The activation energies are 4.66 eV and 5.18 eV 
for oxygen and silicon, respectively (bold solid lines). These numbers compare very well 
with the ones determined in experiments at significantly lower temperatures, namely 4.7 eV 
for oxygen [|D| and 6 eV for silicon [O]. That such an excellent agreement between the the- 
oretical and experimental activation energies is not trivial has recently been demonstrated 
by Hemmati and Angell ||21|| , who showed that various models for silica can give rise to quite 
different activation energies. These authors also showed that the different models predict 
diffusion constants which differ by up to two decades at temperatures as high as 3000 K, 
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which shows that dynamical quantities like D depend much more sensitively on the potential 
than structural quantities. 

For higher temperatures we see from the figure that significant deviations from the Ar- 
rhenius behavior are observed in that the diffusion constants increase slower with increasing 
temperature than expected from an activated process. A similar change in the temperature 
dependence of D has also been found in many other models for silica Recently Hess et 
al. have reported the analysis of experimental viscosity data of SiC>2 at high temperatures 



and find that deviations from a pure Arrhenius law are present ||42|| . 

One explanation to rationalize the observed deviations from the Arrhenius behavior 
is offered by MCT 0, since the so-called ideal version of this theory predicts that the 
temperature dependence of the diffusion constant, as well as the inverse of the a— relaxation 
time t(T), is given by a power-law, i.e. 

D oc t- 1 oc (T - T c y , (6) 

where the critical temperature T c and the critical exponent 7 can, in principle, be calculated 
from the temperature dependence of the partial structure factors. In practice, however, the 
two quantities are usually taken as fit parameters U (for exceptions see j|3| and references 



therein). Instead of using the diffusion data alone to determine 7 and T c , it has been found 
in similar types of analysis (e.g. [|18H20|J44}1 ) that more reliable results are obtained if one 



determines these quantities from simultaneous fits to the diffusion constants as well as to 
the a— relaxation times for different wave-vectors, and therefore we proceeded in this way 



too P5 |. The result of these fits is that the critical temperature is 3330 K. This value of 
T c is in excellent agreement with the one determined by Hess et al. in their analysis of 
viscosity data of real silica, which is 3221 K [fi2]| . (We note that we have learned about the 
results of Hess et al. only after having determined our value of T c ). Using our value of T c , we 
obtain for the critical exponent 7 2.15 and 2.05 for the silicon and oxygen diffusion constant, 
respectively. Note that the theory predicts that the value of 7 should be independent of the 
species and the fact that the two values we find are so close together supports this prediction. 
The so obtained power-law fits are included in Fig. § as well and we see that they give a 
good description of the diffusion constants over about 1.5 decades in D. This range is 
significantly smaller than the one found in the case of simple liquids for which the power- 



law can be observed over about three decades |[17|J18|| . As substantiated below, the reason 
for this relatively small range is very likely the fact that in Si02 the relaxation dynamics at 
low temperatures is dominated by strong hopping processes, i.e. jump like motions of the 
ions. Since the ideal version of MCT, i.e. the version that predicts the power-law given by 
Eq. (H), does not take into account these types of processes, it is thus not surprising that this 
version of the theory is applicable only in a quite limited temperature range. The solution to 
this problem might be the so-called extended MCT, i.e. that version of the theory in which 



such hopping processes are taken into account Since, however, so far the details of 

this theory for the a-relaxation have not been worked out, no test can presently be made. 
This situation is different in the f3— relaxation regime, where some of the predictions of the 
extended MCT have been worked out and as we will demonstrate in a different place [JS|, 
that these predictions hold very well for the present system. 

Although Fig. ^ suggests that MCT might be useful to describe the dynamics of the 
present system, the power-law fits in Fig. 13 are of course not a proof that the extracted 
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parameters T c and 7 have a deep physical meaning. The significance of the theory for 
SiC>2 becomes only obvious if also the a- and /3-relaxation dynamics of the time correlation 
functions are analysed carefully, which is done in Refs. [^,0. We mention already at this 
place, however, that the value for the critical exponent found for the diffusion constant, 
7 m 2.1 is slightly smaller than the one found for the a- relaxation time r for larger wave- 



vectors, which is around 2.35 [f45|. The value of this latter exponent is compatible with the 
relation proposed by MCT between these critical exponents and the value of the so-called 
von Schweidler exponent of the /3-relaxation regime. Hence we conclude that the critical 
exponent of the diffusion constant is smaller than expected from MCT, an observation which 
is in agreement with the one made for simple liquids |IB||i7f . 

Having now presented the temperature dependence of the diffusion constants in the 
temperature range of our simulation, it is instructive to compare these results with the one 
of real experiments. For this we show in Fig. |] the diffusion constants as measured by us 
and the ones determined in experiments by Mikkelsen [fHJ and by Brebec et al. f¥Tf. Also 



included are the Arrhenius fits to our data. We see that the extrapolation of these fits to 
the temperature range which is accessible to the experiments overestimates the diffusion 
constant by about one decade in the case of oxygen and by about two decades in the case 
of silicon. From such a plot we can also estimate the temperature at which the BKS model 
would show the experimental glass transition temperature T g , if we would have access to 
computers which are by about a factor of 10 10 faster. For this we read off the experimental 
values of the diffusion constants at the experimental T s =1450 K (horizontal dashed lines) 
and determine the temperatures at which the Arrhenius extrapolation of our data crosses 
these values. From this construction we obtain a T g ^ sim of 1380 K and 1303 K for oxygen 
and silicon, respectively (vertical dashed lines). Thus we come to the remarkable conclusion 
that using the BKS potential it is possible to estimate the experimental glass transition 
temperature to within 10%. 

Another important transport quantity is the shear viscosity rj. Since rj is a collective 
quantity, it is quite demanding to measure it in a simulation with high accuracy. For the 
case of silica the only simulation we know in which rj has been determined is the one by Barrat 
et al. in which the pressure and temperature dependence of the viscosity was determined at 
relatively high temperatures jP| . 



As in the case of the diffusion constant there are two possibilities to calculate 77: from a 
Green-Kubo relation and a generalized Einstein relation [fi9[| . The Green-Kubo relation is 



V = 7^7 J o dt(A aP (t)A aP (0)) , (7) 



where the off-diagonal elements of the pressure tensor are given by 



N N N 

i Q/3 = E^>f + EE^?S (8) 

i=l i=l j>i 



Here F?j is the a component of the force between ions i and j. 
The Einstein-formula reads: 



— !- -lim {{A ap {t) - A aP {Q)f) (9) 
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with 

N 



A af3 (t) = j2mvm<(t) . (io) 



i=l 



Allen et al. have pointed out that in order to obtain correct results it is important to 
represent the term A aj3 {t) — A a p{$) in Eq. (|9[) in a way which is independent of the coordinate 
system [KOI. One possible choice is 



A a/3 {t) - A a(3 (0) = f dt'A aP {t'). (11) 

Jo 

We also mention that care has also been taken in the evaluation of the forces occurring 
in Eq. (W), since it is not possible to use just the ones obtained from the Ewald sums. More 



details on this problem can be found in Refs. [51]. In the course of our calculation we found 



that the two methods discussed to calculate rj give results with similar accuracy, and thus 
can be considered to be equivalent from a numerical point of view. 

Because of the mentioned collective nature of rj, the usual two runs we did at every tem- 
perature were not sufficient to allow to determine rj with a satisfactory statistical accuracy. 
Therefore we performed for all temperatures above 2900 K 18 extra runs in order to calculate 
r). Because of the long relaxation times at T = 2900 K and T = 2750 K it was not possible 
to do that many runs at these temperatures and thus no results for rj have been obtained. 



In Fig. [10] we show the temperature dependence of the viscosity in an Arrhenius plot 
(filled squares in main figure). We see that, similar to our results for the diffusion constants, 
also this transport quantity shows at low temperatures an Arrhenius behavior which crosses 
over to a weaker temperature dependence with increasing temperature. In the temperature 
regime in which the Arrhenius law is observed the data is described very well by an Arrhenius 
law. If we fit the six lowest temperatures with such a law we find an activation energy of 
5.19 eV, which is in very good agreement with the experimental value of 5.33 eV |34|| . By 
using such an Arrhenius law to extrapolate our data to lower temperatures we see that 
for temperatures around 2800 K such an extrapolation underestimates the viscosity of real 
silica, as measured by Urbain et al. |34| (open circles), by about a factor of 10. Thus we 
find that the BKS model underestimates the viscosity of real silica but that the error lies 
in the prefactor and not in the activation energy. (We also mention that there seems to be 



an uncertainty of the experimental value of the prefactor of about a factor of two [52], thus 
it might be that the actual discrepancy between our simulation and reality is less than just 
stated.) We also mention that if the Arrhenius law found for our low temperature data is 
used for estimating the temperature at which the viscosity of our system reaches the value 
10 13 Poise, i.e. the value of the experimental glass transition temperature, we find that this 
value is at 1310 K. Thus we find, in agreement with our results from the diffusion constant, 
that the BKS model underestimates the experimental glass transition temperature by about 
10%. 

Having an independent measurement of the diffusion constants and the viscosity, it is 
possible to check the validity of the Stokes-Einstein relation in this system. This relation is 
given by 

— ^— = A = const , (12) 

TjD 
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where the constant A has, in the phenomenological Eyring theory, the meaning of a length of 
an elementary diffusion step | 53fl . To check whether the left hand side of Eq. (|i~2"D is indeed 
constant we plot its temperature dependence in the inset of Fig. [H]. From this inset we see 
that A depends in the whole temperature interval on T in that it changes from values around 
20A at high temperatures to values around 5A at the lowest temperatures. This finding is 
not surprising in view of the fact that the activation energy of the viscosity (5.19 eV, see 
Fig. [TUp is very close to the one of the diffusion constant of silicon (5.18 eV, see Fig. ^|), 
thus showing that for silicon the product r/Dsi is constant, whereas in the Stokes-Einstein 
relation an additional factor T is present. Thus we conclude that in the temperature range 
considered here, the Stokes-Einstein relation is not a good way to convert viscosity data into 
diffusivities or vice versa. 

From our data on the temperature dependence of the diffusion constants and the viscosity 
we have evidence that with increasing temperature the relaxation behavior of the system 
crosses over from one similar to an activated process to one which is, potentially, described 
well by the ideal version of MCT. On the other hand, the ideal version of the theory predicts 
a power-law divergence to infinity of D and rj at T c , which in our case is at 3330 K, while 
Figs. |8| and [H] show that this is not the case. Thus the need for the extended version of 
MCT is evident. Further support for a change of the transport mechanisms is obtained 
by investigating the temperature dependence of the life time of a bond between a silicon 
and an oxygen atom. From Fig. we see that even at the highest temperature there is a 
well defined minimum between the first and second neighbor peak in the radial distribution 
function g(r) for the Si-0 correlation. Hence it is quite natural to make the definition that 
a silicon and an oxygen atom are bonded if their separation is less than the location of this 
minimum, which we located at 2.35A. In Fig. 11 we show Ps(t), the probability that a bond 
which was present at time zero is still present at time t. From this figure we see that the 
decay time of this probability increases fast with decreasing temperature. More interesting 
is the observation that at low temperatures -Pb(£) does not decay to zero within the time 
span of our simulation, but that the curves end at around a value of 0.2. In Fig. [7| we have 
shown that our simulations have been long enough to allow to observe the diffusive regime 
in the mean squared displacement and in Refs. |26|]45|] we show that also the structural 
correlation functions, such as the intermediate scattering functions F s (q,t) and F(q,t) at 
the FSDP [[49], decay to zero within the time span of our simulation. Hence we conclude 
that in order to become diffusive, or for the decay of the mentioned structural correlation 
functions, it is not necessary that all of the bonds between silicon and oxygen breakQ. 

From the figure one also sees that the shape of the curves does not seem to depend on 
temperature. This can be checked by plotting Psif) versus t/r B , where the "life time" t b (T) 
is defined by requiring that Pb{t~b) = e _1 . In the inset of Fig. [11] we show this type of 
plot and we see that the curves for the different temperatures fall indeed very well on top 
of each other. At short and intermediate times the master curve is approximated well by 
a Kohlrausch- Williams- Watts curve (KWW), exp(— (t/r B ) /3 ), with j3 = 0.90 (dashed line). 



2 We mention that this is the case at high temperatures also. However, since at these temperatures 
the relaxation times r are short our runs were so long, compared to r, that in the time span of the 
simulation even -Pb(£) decays to zero. 
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decay at short times and a power-law decay at long times. 

In order to find out how the breaking of the bonds is related to the diffusion process 



we plot in Fig. [12] the product of the life time tb and the diffusion constants versus the 
inverse temperature. From this figure we recognize that the product tbDq is essentially 
independent of temperature, demonstrating that the breaking of the bond is indeed related 
to the elementary diffusion step of the oxygen atoms. For the silicon atoms the situation 
is more complicated in that the product TsD§i is constant at high temperatures but then 
starts to decrease at low temperatures. Thus we have evidence that at high temperatures 
the diffusion mechanism for the silicon atoms is very similar to the one of the oxygen atoms 
and is likely governed by the collective motion described by MCT. When the system enters 
the Arrhenius regime the situation changes in that now the silicon atoms see a quite different 
local (mean) potential than the oxygen atoms (since the former sit in deeper minima) and 
thus their diffusion motion is not directly related to the breaking of the bond with one of 
its neighbors. 

A more detailed picture of the change of the transport mechanism can be obtained by 
investigating the self part of the van Hove correlation function G s (r,t) |49| which is defined 
by 

i N a 

G^r,t) = —^(5(r-\r l (t)-r t (0)\)) « G {Si, 0} . (13) 

a i=l 

Thus Airr 2 G"(r, t) is the probability to find a particle at time t a distance r away from the 
place it was at t — 0. In Fig. [13] we show this probability for different times for T = 6100 K 
and T = 2750 K. From the figures at the higher temperature we see that the space and 
time dependence of 47rr 2 G" is very regular in that with increasing time the location of the 
peak moves continuously to larger distances. At this temperature no significant qualitative 
difference between the curves for the silicon and oxygen atoms is observed thus showing that 
the transport mechanism of the two types of particles is very similar. 

At low temperatures the situation is quite different. After the ballistic motion of the 
ions, during which the r and t dependence of is qualitatively similar to the one at high 
temperatures, the distribution function shows a peak whose location depends only weakly 
on time. Such a behavior is well known from studies of the dynamics of supercooled liquids 
and is a manifestation of the cage effect discussed above |T7],T!5]. With increasing time the 



height of this peak decreases and to the left of the peak a long tail is observed in the case 
of silicon and a small peak in the case of oxygen (around 2.6 A). This peak can be seen 
better in a lin-log representation of the curves, which is shown in the insets of the figures. In 



their simulation of a soft-sphere system Roux et al. ]17| related the existence of such a peak 
to the presence of hopping processes and thus we have direct evidence that such processes 
exist in the present system as well. Also included in the figure for the oxygen is the location 
of the first peak in the radial distribution function of the 0-0 correlation at 2.6A (vertical 
line), and we see that this location coincides with the location of the secondary peak in 
47rr 2 G®(r, t). (We also mention that most of the oxygen atoms jumping are defects in that 
they are only one- or three- fold coordinated |30fl .) 

For silicon the situation is different in that no secondary peak is visible at any time. 
However, a careful inspection of the curves reveals that at around 1.6A, marked in the 
inset by a vertical line, a change in the slope of the curves can be noticed. This distance 
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corresponds to the Si-0 bond in one tetrahedron. From Fig. |3]a we see that there are still 
quite a few five-fold coordinated silicon atoms in the melt and we have found that at low 
temperatures the transport of the particles often involves such five-fold coordinated silicon 
atoms [[30|| . Thus a typical bond breaking process involves that a one-fold coordinated 



oxygen atom attaches itself to a four-fold coordinated silicon atom, thus moving the latter 
by a distance of the order of half the Si-0 bond-length. Shortly after this, the now five- 
fold coordinated silicon atom breaks the bond with one of its five nearest neighbors, thus 
creating a new one-fold coordinated oxygen atom, and moves again on the order of half a 
bond-length. Thus in this whole process the silicon atom has moved on the order of the 
length of a Si-0 bond and hence 47rr 2 G^(r, t) shows at this distance the mentioned feature. 

Finally we mention that we have also investigated the space and time dependence of 
the distinct part of the van Hove correlation function, Gd(r,t) p9|j49|] and found that this 
dependence is compatible with our discussion of the self part G% (r, t). This is demonstrated 
in Fig. where we show Ga(r,t) for the 0-0 pair at 6100 K and 2750 K. From Fig. FJa 
we see that even at T = 6100 K Gd{r,t) shows at intermediate times a small peak at the 
origin, r = 0, thus showing that the correlation hole that is observed at time zero is filled 
up not only by a continuous influx from other oxygen atoms but to some extent also by 
particles which hop directly into the position of the oxygen atom which was there at time 



zero. This type of dynamics is not observed in simple liquids at high temperatures [16-18 



and is thus probably attributable to the fact that silica forms a well defined network even at 
high temperatures thus making the correlation hole relatively stable and hence permitting 
that a new oxygen atom can jump into it. We emphasize, however, that at these high 
temperatures this jump mechanism is not the dominant transport mechanism, since the self 
part of the van Hove function does not show a secondary peak (see Fig. |T3|b). 

In Fig. [TJ]b we see that at low temperatures the peak at the origin is very pronounced at 
intermediate times and is even observable at t — 15.4 ns, i.e. at times at which the particles 
show a diffusive behavior (see Fig. 0b). This shows that even on these very long time scales 
the local structure of the network still has some resemblance to the one at time zero and 
that the motion of the atoms is indeed given by the hopping of the atoms in this structure. 



IV. SUMMARY AND CONCLUSIONS 

In this paper we have discussed the results of a large scale molecular dynamics computer 



simulation of a silica melt. The potential used is the one proposed by van Beest et al. [ 2"2( 



which has been found in previous simulations 23] to be quite reliable with respect of repro- 



ducing structural properties of amorphous silica at temperatures below the glass transition 
temperature. In the present work we focus on the structural and dynamical quantities in the 
equilibrium melt. Apart from thermodynamic quantities like the pressure and the compress- 
ibility, we have determined the partial structure factors. From these functions we conclude 
that the melt is substantially ordered even at temperatures as high as 6100 K and shows 
intermediate range order (first sharp diffraction peak) already at around 4000 K. Evidence 
that this high degree of local ordering is not just an artifact of our model is given firstly by 



the fact that at these temperatures the viscosity of real silica is surprisingly high |3"4| , thus 
showing that the motion of the atoms is strongly hindered because of the cage-effect, and 
secondly by our finding that, if we cool the system to a temperature at which the structure 
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can be measured, we find that the structure predicted by the simulation agrees very well 
with the experimental neutron scattering function S n (q), thus showing that the structure of 
the model is very realistic. 

From the mean squared displacement of the particles we calculate the diffusion constants 
D and find that at low temperatures they show an Arrhenius dependence with activation 
energies which agree very well with the ones found in experiments. If this temperature 
dependence is extrapolated to values of D which correspond to the diffusion constants of 
real silica around the experimental glass transition temperature T g , we find that these values 
of D are at temperatures which agree with T g to within 10%. A similar result is obtained 
from the viscosity data. Thus we conclude that the model used is very reliable to predict 
the relaxation behavior of real silica and is even able to predict the real glass transition 
temperature with surprising accuracy. 

For temperatures higher than 3200 K we find significant deviations from the Arrhenius 
dependence of D. In that temperature range the diffusion constants are fitted much better by 
a power- law, D oc (T — T c ) 7 , a temperature dependence which is often found in simple liquids 
and which has been proposed by mode-coupling theory (MCT). The critical temperature 
T c is 3330 K, in excellent agreement with extrapolations by Hess et al. of experimental 
data [12]. A similar conclusion can be drawn from the temperature dependence of the 
viscosity. Thus we find that the relaxation dynamics of this system shows a cross-over from 
a dynamics at high temperatures which can be described by the ideal version of MCT to 
an activated dynamics at low temperatures. The existence of this cross-over can also be 
seen in the temperature dependence of the product of the diffusion constants and the life 
time of a Si-0 bond. This product is constant in the high temperature regime for silicon 
and oxygen, whereas it shows a significant (Arrhenius) temperature dependence for silicon 
at low temperatures, hence giving evidence for a change of the diffusion mechanism with 
decreasing temperature. 

As it has been shown elsewhere [B3|, even the activated dynamics at low temperatures 
can be understood to a good part within the framework of MCT, since predictions of the 
theory like the factorization property, or the existence of the von Schweidler law || seem to 
hold very well. Thus we come to the conclusion that one of the main differences between 
strong and fragile glass formers is that in the former the so-called hopping processes, which 
invalidate some of the predictions of the ideal version of MCT, are, at T c , very pronounced 
in strong glass formers. Due to the strong presence of these hopping processes important 
predictions of the ideal MCT, like the presence of a power-law dependence of the transport 
coefficients, are valid only in a very restricted temperature range. Therefore it is important 
that in the analysis of the relaxation dynamics of strong glass formers those predictions of 
MCT are checked which are not affected by the presence of such hopping processes, if one 
wants to make a real test whether or not MCT is able to describe the relaxation dynamics 
of such a system. 
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FIG. 1. Temperature dependence of the pressure. 
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FIG. 2. Radial distribution functions for different temperatures, a) Si-Si, b) Si-O, c) 0-0. 



19 




2500 3000 3500 4000 4500 

T[K] 




2500 3000 3500 4000 4500 



T[K] 



20 



o 
Q_ 



CO 
CO 

Q_ 



0.5 
0.4 
0.3 
0.2 
0.1 



0.0 
2500 



1.0 
0.8 
0.6 
0.4 
0.2 
0.0 



z=6 



(c) 




z=9 




3000 3500 



4000 4500 
T[K] 





(d) 


z=4 




z=5 




z=3 


— *==^^> 





2500 3000 3500 4000 4500 

T[K] 

FIG. 3. Temperature dependence of the partial coordination numbers z. The dashed lines in 
a) and b) are fits with an Arrhenius law. a) Psi-O) b) Po-Sij c) Pq-O, d) Psi-Si- 
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FIG. 4. Partial structure factors S a p{q) for different temperatures. The curves for 300 K were 
obtained by quenching the system from a high temperature state and are therefore not equilibrium 
curves, a) Si-Si, b) Si-O, c) O-O. 
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FIG. 5. Comparison of the neutron scattering function from our simulation (solid line) with 
the experimental data of Price and Carpenter |36| (circles). 
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FIG. 6. Main figure: Value of the total structure factor at q = as a function of temperature. 
Inset: temperature dependence of the compressibility. 
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FIG. 7. Time dependence of the mean squared displacement for different temperatures, a) 
silicon, b) oxygen. 
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FIG. 8. Arrhenius plot of the diffusion constants. Bold solid lines: Arrhenius fits to the data 
at low temperatures with the stated activation energies. Dashed lines: Results of a fit to the high 
temperature data with the power-law predicted by MCT with a critical temperature of 3330 K. 



25 




1.5 2.5 3.5 4.5 5.5 6.5 7.5 

10 4 /T [HC 1 ] 

FIG. 9. Arrhenius plot of the diffusion constants as determined from our simulation (open 
circles, oxygen; filled squares, silicon) and the experimental values for oxygen (bold line) |^0| and 
silicon [^jjl . Thin solid lines: Extrapolation of our data to low temperatures and extrapolation of 
the experimental data to high temperatures. See text for the discussion of the other lines. 
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FIG. 10. Main figure: Arrhenius plot of the shear viscosity from the simulation (filled squares). 
The dashed line is a fit with an Arrhenius law to our low-temperature data. The open circles are 
experimental data from Urbain et al. Inset: Temperature dependence of the left hand side of 



Eq. (12) to check the validity of the Stokes-Einstein relation. 
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FIG. 11. Time dependence of Pb, the probability that a bond between a silicon and an oxygen 
atom which exists at time zero is also present at time t, for all temperatures investigated (main 
figure). Inset: The same data versus the rescaled time i/re, where tb(T) is the decay time of Pb- 
The dashed line is a fit with a KWW-law with (3 = 0.90. 
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FIG. 12. Temperature dependence of the products of the diffusion constants with the decay 
time of Pb- 
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FIG. 13. Space time dependence of the self part of the van Hove correlation function, a) silicon 
at 6100 K. b) oxygen at 6100 K. c) silicon at 2750 K. The vertical line in the inset corresponds 
to 1.6 A, the length of a Si-0 bond, d) oxygen at 2750 K. The vertical line at 2.6 A in the inset 
corresponds to the nearest neighbor 0-0 distance. In panels a) and b) the times increases by about 
a factor of two from one curve to the other. In panels c) and d) the times are 0.015 ps, 0.033 ps, 
0.38 ps, 174 ps, 362 ps, 590 ps, 817 ps, 1.04 ns, 1.2 ns, 4.5 ns, and 15.4 ns. 
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FIG. 14. Space and time dependence of the distinct part of the van Hove correlation function 
for the 0-0 pair, a) T = 6100 K. b) T = 2750 K. In panel b) the curves that are not labeled 
correspond to times 174 ps, 362 ps, 590 ps, and 1.2 ns. 
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